clear; clc; close all;

k = 3;
h = 0.01;
T = 600;
transient = 550;
N = round(T / h);
N_trans = round(transient / h);

mu_list = linspace(-3*pi, 3*pi, 200);
eta_list = [-2*pi, 0, 2*pi];
colors = [0.7 0 0.7; 0 0 0; 0 0.5 0]; % 紫色，黑色，绿色

figure; hold on; box on;

for idx_eta = 1:length(eta_list)
    eta = eta_list(idx_eta);
    v_traj_all = cell(length(mu_list),1);
    
    % 并行积分，计算后期v轨迹
    parfor i_mu = 1:length(mu_list)
        mu = mu_list(i_mu);
        x = [1e-9; 0; 0; mu; eta];
        f = @(x) system_ode(x,k);
        
        % 丢弃暂态
        for t = 1:N_trans
            x = RK4_step(f, x, h);
        end
        
        v_values = zeros(N - N_trans,1);
        for t = 1:(N - N_trans)
            x = RK4_step(f, x, h);
            v_values(t) = x(5);
        end
        
        v_traj_all{i_mu} = v_values;
    end
    
    % 串行提取峰值并收集
    mu_points = [];
    vpeaks = [];
    for i = 1:length(mu_list)
        mu = mu_list(i);
        pks = findpeaks(v_traj_all{i});
        mu_points = [mu_points; repmat(mu, length(pks), 1)];
        vpeaks = [vpeaks; pks];
    end
    
    % 绘制该eta对应所有峰值散点
    scatter(mu_points/pi, vpeaks/pi, 3, colors(idx_eta,:), 'filled');
end

vpeaks_pi = vpeaks / pi;

scatter(mu_points/pi, vpeaks_pi, 8, colors(idx_eta,:), 'filled');

xlabel('\mu', 'FontSize', 12);
ylabel('v peak', 'FontSize', 12);

xticks(-3:1:3);
xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'});

yticks(-4:1:4);
yticklabels({'-4\pi', '-3\pi', '-2\pi', '-\pi', '0', '\pi', '2\pi', '3\pi', '4\pi'});

grid on;

legend({'\eta = -2\pi', '\eta = 0', '\eta = 2\pi'}, 'Location', 'best');


%% 子函数：系统方程
function dx = system_ode(x, k)
    % 状态变量 x = [x; y; z; u; v]
    dx = zeros(5,1);

    dx(1) = x(2) + x(3) - k * cos(x(5)) * x(2); % \dot{x}
    dx(2) = -x(1) + x(3);                        % \dot{y}
    dx(3) = -x(1) - x(3) + k * cos(x(4)) * x(1); % \dot{z}
    dx(4) = x(1);                               % \dot{u}
    dx(5) = x(2);                               % \dot{v}
end

%% 子函数：RK4一步积分
function x_next = RK4_step(f, x, h)
    k1 = f(x);
    k2 = f(x + h/2 * k1);
    k3 = f(x + h/2 * k2);
    k4 = f(x + h * k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
